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Abstract. In this contribution we consider localized, robust and efficient a-posteriori error 
estimation of the localized reduced basis multi-scale (LRBMS) method for parametric elliptic 
problems with possibly heterogeneous diffusion coefficient. The numerical treatment of such 
parametric multi-scale problems are characterized by a high computational complexity, arising 
from the multi-scale character of the underlying differential equation and the additional param¬ 
eter dependence. The LRBMS method can be seen as a combination of numerical multi-scale 
methods and model reduction using reduced basis (RB) methods to efficiently reduce the com¬ 
putational complexity with respect to the multi-scale as well as the parametric aspect of the 
problem, simultaneously. In contrast to the classical residual based error estimators currently 
used in RB methods, we are considering error estimators that are based on conservative flux 
reconstruction and provide an efficient and rigorous bound on the full error with respect to the 
weak solution. In addition, the resulting error estimator is localized and can thus be used in 
the on-line phase to adaptively enrich the solution space locally where needed. The resulting 
certified LRBMS method with adaptive on-line enrichment thus guarantees the quality of the 
reduced solution during the on-line phase, given any (possibly insufficient) reduced basis that 
was generated during the offline phase. Numerical experiments are given to demonstrate the ap¬ 
plicability of the resulting algorithm with online enrichment to single phase flow in heterogeneous 
media. 
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1 Introduction 


We are interested in efficient and reliable nnmerical approximations of elliptic parametric 
mnlti-scale problems. Snch problems consist of finding p(/i) G Q, snch that 

y) = Kq) for all qeQ, (1) 

in a snitable space Q for parameters fx € V G for p G N, where the data fnnctions 
involved may depend on an a-priori given mnlti-scale parameter e > 0 (described in detail 
in Sectionj^. An approximation Phipi-) G Qh of p{fx) is nsnally obtained by a discretiza¬ 
tion of Q resnlting from a Galerkin projection onto a discrete space Qh, associated with 
a triangnlation Th of Since Th shonld resolve all featnres of Q associated with the 
fine scale e, solving parametric heterogeneons mnlti-scale problems accnrately can be 
challenging and compntationally costly, in particnlar for strongly varying scales and pa¬ 
rameter ranges (see the references in j44jl. Two traditional approaches exist to rednce the 
compntational complexity of the discrete problem: nnmerical mnlti-scale methods and 
model order rednction techniqnes. Nnmerical mnlti-scale methods rednce the complex¬ 
ity of mnlti-scale problems with respect to e for a fixed /x, while model order rednction 
techniqnes rednce the complexity of parametric problems with respect to /x for moderate 
scales e (see m for an overview). In general, nnmerical mnlti-scale methods captnre the 
macroscopic behavior of the solntion in a coarse approximation space Qh C Qh, e.g., 
associated with a coarse triangnlation Th of fl, and recover the microscopic behavior of 
the solntion by local fine-scale corrections. Inserting this additive decomposition into Q 
yields a conpled system of a fine- and a coarse-scale variational problem. By appropriately 
selecting trial and test spaces and defining the localization operators to deconple this sys¬ 
tem, a variety of nnmerical mnlti-scale methods can be recovered, e.g., the mnlti-scale 
finite element method (MsFEM) |22) . the variational mnlti-scale method |43| 137] . the 
mnlti-scale finite volnme method |30| and the heterogeneons mnlti-scale method (HMM) 
13 EU, jnst to name a few. Model order rednction nsing rednced basis (RB) methods, on 
the other hand, is based on the idea to introdnce a rednced space Qred C Qh, spanned 
by discrete solntions for a limited nnmber of parameters /x. These training parameters 
are iteratively selected by an adaptive Greedy procednre (see e.g. |56j and the references 
therein). Depending on the choice of the training parameters and the natnre of the prob¬ 
lem, Qred is expected to be of a significantly smaller dimension than Qh- Additionally, if 
b allows for an affine decomposition with respect to /x, its components can be projected 
onto Qred) which can then be nsed to effectively split the compntation into an off-line and 
on-line part. In the off-line phase all parameter-independent qnantities are precompnted, 
snch that the on-line phase’s complexity only depends on the dimension of Qred- The 
idea of the recently presented localized rednced basis mnlti-scale (LRBMS) approach 
(see mm is to combine nnmerical mnlti-scale and RB methods and to generate a local 
rednced space C for each coarse element of Th, given a tensor prodnct type 
decomposition of the fine approximation space, Qh = ©tsThQ^- The coarse rednced 
space is then given as QrediTn) '■= ^TeTnQjed Qh, resnlting in a mnltiplicative de¬ 
composition of the solntion into Pred{x; /x) = Pnix] pL)ipn{x), where the RB 
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functions (pn capture the microscopic behavior of the solution associated with the fine 
scale £ and the coefficient functions pn only vary on the coarse partition Th- 

Other model reduction approaches for parametric (but not multi-scale) problems that 
incorporate localization ideas and concepts from domain decomposition (DD) methods 
are the reduced basis element method |391 [38] , the reduced basis hybrid method |331 IM] 
and the port reduced static condensation reduced basis element method |52j . While the 
idea of the former is to share one reduced basis on all subdomains the idea of the latter 
two is to generate one reduced basis for each class of subdomains which are then coupled 
appropriately. There also exist several approaches to use model reduction techniques 
for homogenization problems (see IH) and problems with multiple scales, such as the 
reduced basis finite element HMM Em- For the case of no scale-separation there exist 
the generalized MsFEM method |21j . which incorporates model reduction ideas, and 
most recently a work combining the reduced basis method with localized orthogonal 
decomposition (see El). 

It is vital for an efficient and reliable use of RB as well as LRBMS methods to have 
access to an estimate on the model reduction error. Such an estimate is used to drive 
the adaptive Greedy basis generation during the off-line phase of the computation and 
to ensure the quality of the reduced solution during the on-line phase. It is usually given 
by a residual based estimator involving the stability constant and the residual in a dual 
norm. It was shown in [8] that such an estimator can be successfully applied in the 
context of the LRBMS, but it was also pointed out that an estimator relying on global 
information might not be computationally feasible since too much work is required in 
the off-line part of the computation. 

The novelty of this contribution lies in a completely different approach to error esti¬ 
mation - at least in the context of RB methods. We make use of the ansatz of local error 
estimation presented in |26] which measures the error by a conforming reconstruction 
of the physical quantities involved, specifically the diffusive flux This 

kind of local error estimation was proven to be very successful in the context of multi¬ 
scale problems and very robust with respect to e. We show in this work how we can 
transfer those ideas to the framework of the LRBMS to obtain an estimate of the error 
|||p(/x) — Pred(/.t)|||^. We would like to point out that we are able to estimate the error 
against the weak solution p(/x) in a parameter dependent energy norm while traditional 
RB approaches only allow to estimate the model reduction error in a parameter indepen¬ 
dent norm and only against the discrete solution. In principal, this approach is able to 
turn the LRBMS method into a full multi-scale approximation scheme, while traditional 
RB methods can only be seen as a model reduction technique. We would also like to 
point out that, to the best of our knowledge, this approach (first published in |45) 1 is the 
first one to make use of local error information in the context of RB methods. 

2 Problem formulation and discretization 

We consider linear elliptic problems in a bounded connected domain 12 C for d = 2,3, 
with polygonal boundary dO, and a multiplicative splitting of the influences of the multi¬ 
scale parameter e and the parameter pL. An example is given by the problem of finding 
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( 2 ) 


a global pressure p(/i) G Hq{Q) for a set of admissible parameters £ V, such that 


-V- {X{fi)KeVp{n)) = f 


in n 


in a weak sense in i?Q(n), where denotes the usual Sobolev space of weakly differ¬ 
entiable functions and Hq{Q) C its elements which vanish on the boundary in 

the sense of traces. Problems of this kind arise for instance in the context of the global 
pressure formulation of two-phase flow in porous media. Using an IMPES discretiza¬ 
tion scheme of the pressure/saturation system (see m and the references therein), an 
equation of the form (|^ has to be solved in each time step for a different parameter 
H £ V. In that context A(/x) denotes the scalar total mobility and denotes a possibly 
complex heterogeneous permeability tensor; external forces are collected in / and the 


parameter /x models the influence of the global saturation (see Definition 2.1 for details 
on A, Kg and /). Note that more complex boundary conditions and additional parameter 
dependencies of the boundary values as well as the right hand side (modeling parameter 
dependent wells, for instance) are possible, but do not lie within the scope of this work. 

Triangulation. We require two nested partitions of D, a coarse one, Th, and a fine 
one, Th- Let Th be a simplicial triangulation of D with elements t £ Th- In the context 
of multi-scale problems we call Th a fine triangulation if it resolves all features of the 
quantities involved in ([^, specifically if £ [L°°(t)]'^^'^ is constant for all t £ Th- 
This assumption is a natural one in the context of two-phase flow problems where the 
permeability field is usually given by piecewise constant measurement data. For simplicity 
we require Th to fulfill the requirements stated in |26l Sect. 2.1], namely shape-regularity 
and the absence of hanging nodes; an extension to more general triangulations is possible 
analogously to |26[ A.l]. We only require the coarse elements T £ Th to be shaped such 
that a local Poincare inequality in H^{T) is fulfilled (see the requirements of Theorem 
4.2). We collect all fine faces in all coarse faces in Th and denote by M{t) C Th and 


W(r) C Th the neighbors oit £ Th and T £ Th-, respectively and by /i* the diameter of 
any element * of the sets Th, Th, Th or Th- We collect in t"^ C Th the fine elements of Th 
that cover the coarse element T and in Th C Th all faces that cover the set *, e.g. by Tj^ 
the faces of a fine element t £ Th,hy Tj^ the faces that cover a coarse face E £ Th and so 
forth; the same notation is used for coarse faces Tfi C Th- In addition we denote the set 
of all boundary faces hy Th C. Th and the set of all inner faces, that share two elements, 

o _ o _ o 

by Th C Th, such that Th D Th = Th and Th D Th = ^- We also denote the set of fine 

_ rp 

faces which lie on the boundary of any coarse element T G Th by Th '-= U_Ee.F^ 
by T'J^ := T^\tJ^ the set of fine faces which lie in the interior of the coarse element. 
Finally, we assign a unit normal Ug G to each inner face dt~ n dt~^ = e £ Th, pointing 
from t~ to t^, and also denote the unit outward normal to dCl by rig for a boundary face 
e = dt~ n 5D, for t~,t~^ £ Th- 


2-1 The continuous problem 

For our analysis we define the broken Sobolev space H^{Th) C LP‘{VL) by H^{Th) := 
{q £ I q\t £ H^{t) Vt G r/i} as the most general space for the weak, the 
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discrete and the reduced solution defined below. In the same manner we denote the 
local broken Sobolev spaces C L‘^{T) for all coarse elements T G Th and denote 

by ; H^{Th) —t [L^(n)]‘^ the broken gradient operator which is locally defined by 
{'Vhq)\f. := V(gU for all t £ and q G Now given A(/x), and / as stated 


in Definition 2.1 we define the parametric bilinear form b : V ^ [H^{Th) x H^{Th) 

M], /I I—)• [{p,q) I—t 6(p, g;/x)], and the linear form I : H^{Th) —M by b{p,q]fi) := 
YIiT&t Kq) StsT respectively, and their local counterparts 

b'^ : V ^ X — )> M], /x i-)- [{p, q) b'^{p, q; /x)], and —)> M for 

all T G Th, £T and p,q G H^{Th) by 


b^{p,q]fj.) := j (A(/x)K£Vfep)-Vfegdx 


and 




/^dx, 


respectively. 

Definition 2.1 (Weak solution). Let f G Lp‘{Vl) be bounded, X{n) G L°°{Ll) be strictly 
positive for all ^ £ V and G symmetric and uniformly positive definite, 

such that A(/i)k£ G is bounded from below (away from 0) and above for all 

/jl £ V. We define the weak solution p :V ^ Hq{Q) for a parameter pi £ V, such that 


Kp{i^),T,tj) = Kq) 


for all q £ Hq{LI). 


(3) 


Note that, since 6(-, •; /x) is continuous and coercive for all /x G "P (due to the assump¬ 
tions on X{pL)Ke) and since I is bounded, there exists a unique solution of ^ due to the 
Lax-Milgram Theorem. Given these requirements we denote by 0 < c],{pL) and < 

the smallest and largest eigenvalue of (A(/x)k£)|j, respectively, for any £ V and 
t £ Th and additionally define 0 < c* ;= min^g-p Cg(/x) and c* < C(, := max^gp C({pL). 

A note on parameters. In addition to the assumptions we posed on A above 
we also demand it to be affinely decomposable with respect to pi £ V, i.e there ex¬ 
ist .S’ G N strictly positive coefficients 0^ : P —?■ M and S nonparametric components 
X^ £ L°°{'Ll), for 1 < ^ < S', such that X{x] ph) = We can then com¬ 

pare A for two parameters /x,/! G P by a{pi,~pL)X{ji) < A(/x) < 7 (/x, 7l)A(7x), where 
a{pi,Ji) := min|L^ 0^(/x)0^(7J)~^ and 'y{pi,Ji) := max|L^ 0^(/x)6l^(7l)~^ denote the posi¬ 
tive equivalence constants. This assumption on the data function A is a common one 
in the context of RB methods and covers a wide range of physical problems. If A does 
not exhibit such a decomposition one can replace A by an arbitrary close approximation 
using Empirical Interpolation techniques (see |in| I2n| ) which does not impact our anal¬ 
ysis. All quantities that linearly depend on A inherit the above affine decomposition in 
a straightforward way. In particular there exist component bilinear forms such that b 
and 6^ (and their discrete counterparts introduced further below) are also affinely de¬ 
composable. Since we would like to estimate the error in a problem dependent norm 
we introduce a parameter dependent energy (semi) norm Ill-Ill. : P —)• [H^{Th) —>■ M], 
P ^ [9 ^ llklll;x]> defined by |||g|||^ := [Y.t&Th IH^IHut)^ i^h the local semi-norm 
defined by IUq'III^'j- := b'^{q,q‘, for all T £ Th, p G P and q £ H^{Th)', the local 
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semi-norm on the fine triangulation, for any t G Th, is defined analogously to 

IIMII/xT- Note that Ill-Ill. is a norm only on We can compare these norms for any 

two parameters € V using the above decomposition of A (the same holds true for 
the local semi-norms): 

\fa{p7p) IIMIIp: < IIMII^ < IIMIIp: (4) 

2.2 The discretization 

We discretize ^ by allowing for a suitable discretization of at least first order inside 
each coarse element T G Th and by coupling those with a Symmetric weighted Interior- 
Penalty discontinuous Galerkin (SWIPDG) discretization along the coarse faces of Th- 
This ansatz can be either interpreted as an extension of the SWIPDG discretization 
on the coarse partition Th, where we further refine each coarse element and introduce 
an additional local discretization, or it can be interpreted as a domain-decomposition 
approach, where we use local discretizations, defined on subdomains given by the coarse 
partition, which are then coupled by the SWIPDG fluxes. In view of the latter, this ansatz 
shows some similarities to [16] but allows for a wider range of local discretizations. A 
similar ansatz for a multi-numerics discretization using a different coupling strategy was 
independently developed and recently introduced in [3^. We will present two particular 
choices for the local discretizations and continue with the definition of the overall DG 
discretization. 

Local discretizations. The main idea is to approximate the local bilinear forms 6^, 
which are defined on the local subdomain triangulations , by local discrete bilinear 
forms 5^ discretizing ([^ on T with homogeneous Neunrann boundary values. We will 
additionally choose local discrete ansatz spaces C with polynomial order 

/c G N, to complete the definition of the local discretizations. The first natural choice 
for the local discretization is to use a standard continuous Galerkin (CG) discretization, 
which we obtain by setting = b'^ and Qj^’ = with 

Sf^iT^):={q€C\T)\q\,eFk{t) Vt G C H\T), 

where lPfc(a;) denotes the set of polynomials on a; C D with total degree at most A: G N. 
Another choice is to use a discontinuous space given by 

Q'^^{Tl):={qeL\T)\q\,eF,{t) Vf G C 

and to choose 6^ from a family of DG discretizations. Therefore we introduce the techni¬ 
calities needed to state a common framework for the non symmetric, the incomplete, the 
symmetric and the symmetric weighted interior-penalty (IP) DG discretization (hence¬ 
forth denoted by NIPDG, IIPDG, SIPDG and SWIPDG, respectively, see [23 and the 
references therein), following |26l Sect. 2.3]. 

For a function q G H^{Th) that is double-valued on interior faces we denote its jump 
on an inner face e £ Thhy [[g]]g := q~ —q'^ with q^ := q\^±, recalling that e = t~ Dt^ for 
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^ Th- We also assign weights uj~, coJ~ > 0 to each inner face, such that uj~ + = 1, 

and denote the weighted average of q by {{^jjg := ^ boundary face 

e G Fh we set uj~ = 1, a;+ = 0, [[g]]g := q and {{gjjg := q. With these definitions we define 
the local bilinear form 6^ : iP — )• x — )■ M], /x [(p, q) bj^{p, q; /x)] for 

r?G {-1,0,1} by 

bl{p,q;f^) ■=b'^{p,q;n)+ /x) + 6^(p, g;/x) + 6^(g,p;/x)) (5) 

on T G Th, with its coupling and penalty parts b^ and 6®, respectively, defined by 

:= y-{{('^(p)'^e'^p)-’^e}}g Mgds and b^{p,q-, fx) := j '7g(/x) [[p]]g [[gj^ ds , 
e e 

for all /X G "P, all p,q G H^{Th) and all e G Fh- The parametric positive penalty function 
crg(/x) : P* —7- M is given by CTeifx) := {{A(/x)}}g (t|, where <t > 1 denotes a user- 

dependent parameter and the locally adaptive weight is given by cj| := 6^ 6~{6^ + 

o _ 

for an interior face e G Fh and by it| := 5~ on a boundary face e G Fh, respectively, 
with 5^ := ne-K^ -ng. Using the weights uj^ := 1/2 we obtain the NIPDG bilinear 
form for ?? = —1, the IIPDG bilinear form for ?? = 0 and the SIPDG bilinear form for 
'!? = !. We obtain the SWIPDG bilinear form for = 1 by using locally adaptive weights 
io~ := + 6~)~^ and ujf := S~(6^ -|- From now on we assume that bj^ is 

of the form ([^, since this is the most general case and also covers a GG discretization, 
where all face terms vanish due to the nature of 5^(r^). 

Global coupling. Now given suitable local discretizations on the coarse elements T G 
Fh we couple those along the coarse faces E G Fh using a SWIPDG discretization again 
and define the bilinear form bh ■ V ^ [H^{Th) x H^{Th) —t M], /x i—)• [{p, q) i—)• bh{p, q] p)], 
by 

bh{p,q]tJ-) ■= Y ^ft(pw;p)+ Y (6) 

T&Th E&Hh 

where we use the SWIPDG variants of to define the coupling bilinear form b^ :V ^ 
[H^{Th) X H^{Th) M], p [(p,g) 6f (p,g;p)], by 

bhipF-,tJ-) ■■= Y (^c(9>p;p)+ ^c(p7 9;p)+ &p(9,p;p)) (7) 

for all E G Fh, all p G P" and all p,q G H^{Th)- Accordingly, we define the DG space 
Q\ C H^{Th) for A: > 1 by 

Qi-.= {qGH\Th) I qW^QY ^TgTh}, 
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with either being the local CG space or the local DG space 

Definition 2.2 (DG solution). Let a be large enough, sueh that bh{-, S/^) is continuous 
and coercive for all £ V. We define the DG solution : V ^ Qj^ for a parameter 
tieV, such that 

hh{ph{ti), Qh', m) = l{Qh) for all qh G Q\. (8) 

As always with DG methods, coercivity is to be understood with respect to a DG norm 
on Qjj (i.e., given by the semi norm combined with a DG jump norm) and there exists 
a unique solution of (|^ due to the Lax-Milgram Theorem, if a is chosen accordingly. 

Remark 2.3 (Properties of the discretization). Depending on the choice of the coarse 
partition and the local discretization we can recover several discretizations from ([^. 
Choosing Th = ^ ond yields a standard CG discretization (except for 

the treatment of the boundary values), for instance. Choosing Th = Th, on the other 
hand, results in the standard SWIPDG discretization proposed in j2'7\ \2(^ . Note that the 
local discretizations as well as the polynomial degree k do not have to coincide on each 
coarse element (or even inside a coarse element when using a local DG space). It is 
thus possible to balance the computational effort by choosing local CG or k-adaptive DG 
discretizations. This puts our discretization close to the multi-numerics discretization 
proposed in where the latter allows for an even wider range of local discretizations 
while coupling along the coarse faces using Mortar methods. Concerning the choice of 
the user dependent penalty factor, we found an automated choice of a depending on the 
polynomial degree k, as proposed in JMEj, to work very well. 

3 Model reduction 

Disregarding the multi-scale parameter e for the moment, model reduction using reduced 
basis (RB) methods is a well established technique to reduce the computational complex¬ 
ity of ([^ with respect to the parameter /x. It has been successfully applied in multi-query 
contexts, where ([^ has to be solved for many parameters /x (e.g. in the context of op¬ 
timization), and in real-time contexts, where a solution of ([^ (or a quantity depending 
on it) has to be available for some fi as fast as possible (see j8l[20ll33] and the references 
therein). The general idea of RB methods is to span a reduced space Qred C Q\ by post- 
processed solutions of (|^ for a selection of parameters and to apply a Galerkin projection 
of ([^ onto this reduced space. Together with the assumption we posed on the data func¬ 
tion in section 1^ this allows for the well-known off-line/on-line decomposition of the RB 
method, where all quantities depending on (in particular high-dimensional evalua¬ 
tions of the bilinear form) can be precomputed and stored in so called Gramian-matrices 
during a possibly computationally expensive off-line phase. In the on-line phase of the 
simulation only low dimensional quantities depending on the parameter dependency of 
the data functions and the dimension of Qred need to be computed which can usually be 
done in real time and even on low end devices such as smartphones or embedded devices. 
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It has been shown in j^, however, that the computational complexity of the off-line 
phase can become unbearable if the RB method is applied to parametric multi-scale 
problems, such as (|^ for small scales e <C 1. There are two main drawbacks of classical 
RB methods in this context: (i) all high-dimensional quantities (solution snapshots, 
operator matrices, ...) usually depend on the fine triangulation Th, the size of which 
scales with 0{l/£). Correspondingly all evaluations involving these quantities (scalar 
products, operator inversions, ...) become increasingly expensive for small scales e; (m) 
the generation of the reduced space Qred usually requires the discrete problem Q to be 
solved for several parameters /x. In the context of multi-scale problems, however, it is 
usually only feasible to carry out very few global solutions of the original problem on the 
fine triangulation, if at all. 

The first of these shortcomings has been successfully addressed by the localized reduced 
basis multi-scale (LRBMS) method introduced in |36[ [8], which takes advantage of the 
coarse partition and the underlying discretization to form local reduced spaces C 
on each coarse element instead of the usual global RB space. This allows to balance 
the computational effort of the off-line versus the on-line phase by selecting an appropriate 
size of the coarse partition Th depending on the multi-scale parameter e and the real-time 
or multi-query context of the application. But it was also shown in [8], that the standard 
residual based estimator which is usually used in the context of RB methods does not 
scale well in the context of multi-scale problems. We thus finalize the definition of the 
LRBMS method in this contribution with an efficiently off-line/on-line decomposable 
a-posteriori error estimator (see the next section). The second shortcoming of classical 
RB methods is addressed in section where we propose an on-line enrichment extension 
for the LRBMS method based on the presented error estimator, which finally turns the 
LRBMS method into a fully capable multi-scale method. 

The reduced problem. Let us assume that we are given a local reduced space 
C on each coarse element T G Th- Those are usually, but not necessarily, 

spanned by solutions of Q and of low dimension, i.e., dimQ^^ <C dimQ^’^. We then 
define the coarse reduced space QrediTn) C Q\ by 

QrediTn) := {q^Qi \ q\T ^ Qld VT G Th] C QI 

and obtain the reduced solution by a Galerkin projection of ([^ onto QrediTn)- 

Definition 3.1 (Reduced solution). We define the reduced solution pred : T —)• QrediTn) 
for a parameter p G "P, such that 

bhiPredifJ-), qred] M) = Kqred) for all Qred G QrediTn)- (9) 

Off-line/on-line decomposition. The well-known off-line/on-line decomposition of 
RB methods relies on an affine parameter dependence of the data functions (see section]^ 
that carries over to the discrete bilinear form, i.e., there exist S nonparametric component 
bilinear forms : H^irQ x H^irQ —>■ M, for 1 < ^ < S’, such that bhip, q', fJ-) = 
Yl^=i^iif-'')bh,^ip,q), with 6^ being the coefficients of the corresponding decomposition 
of A. In the standard RB framework, given a basis {qred,i, - - -, qred,N} of Qredj with 
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N := dimQred) one would precompute dense Gramian matrices bh,^ G K"'"'"', given 
by {bh,^)^j •= bh,^{qred,i, Qred,j) during the off-line phase. In the on-line phase, given a 
parameter G V, one would obtain the reduced system matrix by a low dimensional 
summation bh{fi) ■= ^ which only involves evaluations of the 

scalar coefficient functions 6^ for the current parameter. The resulting dense reduced 
system is of size N x N and does not depend on the dimension of Q^. It can thus be 
solved with a complexity of 0{N^) in the on-line phase. Denoting by '■= dimQ^ and 
:= dim(5|(r^) (with = Y^t&Th ^h) number of Degrees of Freedom of the 
global and the local spaces, respectively, the reduction of bh during the off-line phase, 
however, would be of a computational complexity of 0{SN‘^Nh), which can become too 
expensive in the context of small scales e. 

The flexible framework of the LRBMS method, however, allows for a more local ap¬ 
proach. Since the affine decomposition of the data function A also carries over to the local 
bilinear forms 6^ and the coupling bilinear forms b^, we define local Gramian matrices 
bl^ G on each coarse element T G Th, given by 

9red,j) T ^ ^ ^/i,^(9red,i) 9red,j)) 




where i, • • •, qj^^ tvt} denotes a basis of with := dimQred- addition, we 
define coupling Gramian matrices bj^'^ G for all coarse elements T G 7h and 

their neighbors S G AA(T) by 


T,S 






bh,^ (^redji) QiedJ ) ■ 


In the same manner, we define local vectors G , given by {F)j := I’^iq^ed^)^ for all 
T G Th- The global Gramian matrices bh,^ G R^^^ and the global vector I G R^, with 
^ ■= YjTGTh are then given by a standard DG mapping (with respect to the coarse 
triangulation) of their local counterparts. The complexity to directly solve the reduced 
system in the on-line phase, d((X]'rg 7 ^ -^^)^)) is by definition larger than for standard 
RB methods. The reduction of bh, however, can be carried out with a complexity of 
roughly 0{S N'^'^N'^), which can be computed significantly faster than in the 

standard RB setting described above (since S\Th\ local computations can be carried out 
in parallel), depending on the choice of Th- 

Remark 3.2 (Properties of the reduced system). The reduced system matrix of the 
LRBMS method is obviously larger than the one we would obtain for standard RB methods 
(given that it scales with \Th\)- On the other hand we obtain a sparse matrix of dense 
blocks (stemming from the local and coupling blocks and the DG mapping) that is usually 
much smaller than the system matrix of the high-dimensional problem. Thus, if \Th\ is 
large, the reduced system can also he solved using sparse iterative solvers with a complexity 

of (-){{Yt£Th 
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4 Error analysis 


Our error analysis is a generalization of the ansatz presented in |26] to provide an esti¬ 
mator for our DG approximation (|^ as well as for our LRBMS approximation (|^. The 
main idea of the error estimator presented in |57[ |26] is to observe that the approximate 
DG diffusive flux —A(/r)KeV/iP/i(/x) is nonconforming while its exact counterpart belongs 
to -ffdiv(f^) C which denotes the space of vector valued functions the diver¬ 

gence of which lies in L^(D). The idea of |57[ [26] is to reconstruct the discrete diffusive 
flux in a conforming Raviart-Thomas-Nedelec space C and compare it 

to the nonconforming one. Their error analysis relies on a local conservation property 
of the reconstructed flux on the hne triangulation Th to prove estimates local to the hne 
triangulation. 


We transfer this concept to the DG discretization dehned in section 2.2 and prove 
estimates local to the coarse triangulation that are valid for the DG approximation 
as well as for our LRBMS approximation. We obtain mild requirements for the coarse 
triangulation and the local approximation spaces, namely that a local Poincare inequality 
holds on each coarse element and that the constant function 1 is present in the local 
approximation spaces. The latter is obvious for traditional discretizations and can be 
easily achieved for the LRBMS approximation. The estimates are fully off-line/on-line 
decomposable and can thus be efficiently used for model reduction in the context of the 
LRBMS. Preliminary results were published in |45j . 

We begin by stating an abstract energy norm estimate (see |26l Lemma 4.1]) that splits 
the difference between the weak solution p G iLg(D) of ^ and any function q G 
into two contributions. This abstract estimate does not depend on the discretization and 
thus leaves the choice of s and v open. Note that we formulate the following Lemma 
with different parameters for the energy norm and the weak solution. The price we have 
to pay for this flexibility is the additional constants involving and that 

vanish if Ji and /x coincide. In the following we denote the product over a space R(w), 
for a; C D, by (•, ■)■[/,oj and omit w if a; = D; the same holds for the induced norm IHly^. 

Lemma 4.1 (Abstract energy norm estimate). Let p{pi) G 77g(D) be the weak solution 
of @ for /X G "P and let ph G H^{Th) and pi be arbitrary. Then 


lb(p) -Phlllu < 


< 

^ _( 

inf 

n — 


seH^in) 

+ 

inf 

sup 


V^ti din 





< 

y/iinT) 

y/oiinp) 

2 IIIp(m) 


( 10 ) 


Proof. We mainly follow the proof of |261 Lemma 4.1] while accounting for the parameter 
dependency of the energy norm and the weak solution. It holds for arbitrary /x G P, 
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p G and ph G H^{Th), that 


■Ph\L< inf |bh-s|| 


+ sup b{p - Ph,(p, ti) 
iiv’iii, =1 


(11) 


(see m Lemma 7.1]) and for the weak solution p{pi) G Hq{^) of (|^, that 
b{p{n) -ph,(f,fJ,) = (/,<^)^2 - (A(Ai)K£Vfep/i, Vv7)^2, 

= (/ - V-U,(/?)^2 - {X{ti)Ks-VhPh + V,Vip) 


( 12 ) 


for all (p G L7g(n) and all v G L7div(fl)) where we used the definition of b in the first 
equality and the fact that {v,Vp )]^2 = —{'V-v,(p)i 2 due to Green’s Theorem and p G 
iLg(fl) in the second one. Inserting (12) into 0 L with p = p{pi) and using the norm 
equivalence (Q then yields the first inequality in (10). 

To obtain the second inequality we choose s = p{fJ^) and v = —X{n)KeVp{fi) in the right 


hand side of (10) which eliminates the two infimums and leaves us with two terms yet 


to be estimated arising inside the supremum. Using Green’s Theorem and the definition 
of b we observe the vanishing of the first term. We estimate the second term as 


\{X{fj,)KeVhPh - A(/x)k£Vp(/x), V(/?)^2| 

= I {{X{tJ-)Key^‘^Vh{ph - p{fi)), (A(/i)Ke)^/^V¥5)^2 

< |(A.... 

-P(/^)III;X 


(A(/x)K£)^/^Vh(ph -P(/^)) {X{tJ,)Ke)^^^Vip 

L‘^ 


L2 


where we used the Gauchy-Schwarz inequality and the definition of the energy norm. We 
finally obtain the second inequality of (10) from the bound above by observing that the 
supremum vanishes (due to IHy^lH^ = 1) and by using the norm equivalence (|^ again. □ 

The following theorem states the main localization result and gives an indication on 
how to proceed with the choice of u: it allows us to localize the estimate of the above 
Lemma, if v fulfills a local conservation property. It is still an abstract estimate in the 
sense that it does not use any information of the discretization and does not yet fully 
prescribe s and v. 

Theorem 4.2 (Locally computable abstract energy norm estimate). Let p{fi) G Hq(Q) 
be the weak solution of (|^ for /x G "P, let s ^ Hq{Q) and ph G H^{Th) be arbitrary, let 
V G fulfill the local conservation property (V -v, 1)^2 = (/, 1)^,2 j’ and let Cp > 

0 denote the constant from the Poincare inequality ||(/9 — ng(/9||^2 ^ < Cp/i^H V (^||^2 p for 
all p G H^{T) on all T G Th, where Hf denotes the Lf -orthogonal projection onto IP/(w) 
for I G N and w C Q. It then holds for arbitrary fl, fi € V, that 


lb(p) - PklW-n; < PiPh, S, V, /X, fl), 
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with the abstract global estimator i){ph, s,v, fj,, fi) defined as 


viPh,s,v;n,fi) := 




Yj vLiPh,s-,t^fj + (^ Y 

T&Th T&Th 


1/2 / ^ 

vY 


+ 75b)( S ’'I*'- 

^ T&Th 


1/2 


and the local nonconformity estimator defined as ri'^fiph, s;Jl) := \\\ph — local 

residnal estimator defined as rjJ(v) := {Cp/c^Y^‘^hT\\f — 'V-v\\p 2 rp and the local diffnsive 
flnx estimator defined as fj^fiph^v^fi) := {^{pY^Y^hPh + v)^p 2 ^ for all 

coarse elements T G Th, where cf := 


miiiterj4- 


Proof. We loosely follow the proof of |261 Theorem 3.1] while acconnting for the parameter 
dependency and the coarse triangnlation. Fixing an arbitrary s G Hq{Q,) in (10) and 
localizing with respect to the coarse triangnlation yields 


Mp) -Phlh < 


.k-'iJY. 


— s 


JI,T 


T&Th 


+ snp { E (/ - V-n,99) 1^2 T {Kp) ■ ^hPh + V,Vip) p2 q 

ip£Hi(U) 1 T'c- 7 -„ ---^ '--- 


(13) 


|||¥>|||^ = 1 :=(*) :={ii) 

which leaves ns with two local terms we will estimate separately. 

(i) Since (/ — V-f, n ^(^)^2 = 0 dne to the local conservation property of v we can 
estimate the first term as 

\{f -V-v,(p)p 2 ^t\ - Wf 

<^f^hT(max^] ||/-V-n||i2; 

Cg/ 




where we nsed the Canchy-Schwarz ineqnality, the Poincare ineqnality and the local 


norm eqnivalence (Q on all t ^ Tp . 
(a) We estimate the second term as 


\{X{fJ.)KeVhPh + V,^p)L2^T\ - iKp)l^e) {Kp)l^e^hPh + v) 
<YYp^ ^ {\{ij)Ke)~^^‘^{\{pi)KeVhPh + v) 


L^,T 


P,T 




M,r 


nsing the Canchy-Schwarz ineqnality, the definition of the local energy semi-norms 
and the parameter eqnivalency from section 


13 

























Inserting the last two ineqnalities in (13) and nsing the Canchy-Schwarz inequality yields 


lb(M) -Ph\\\T[< 






T&Th 


1/2 


7ii'o\ L j~:z 


+ sup _ 

ip£H^{n) ^TgTh 

Ill9=lt=l ^ 


v) + 


\/a(M.A) 


Vdf{Ph,V,(l 


(Hi) 


1 2 


/x,T 


using the definition of the local estimators and of . Using the Canchy-Schwarz inequal¬ 
ity again we can further estimate (in) as 


(in) < 


E 

TgTh 


~T 

Vr 


{v)‘ 


1/2 


+ 


\/a(M.A) 


E 

T£Th 


vli{ph,v,fif 


1/2 


The previous two inequalities combined give the final result, since the supremum vanishes 


due to 


= 1 . 


□ 


Remark 4.3 (Properties of the locally computable abstract energy norm estimate). In 
contrast to the estimator proposed in the above estimate is local with respect to Th, 
not T/j. Choosing 7 h = Th we obtain nearly the same estimate as the one in 126^ for the 
pure diffusion case (apart from a slightly less favorable summation). In general, however, 
we can only expect fj,- to be superconvergent if we refine Th along with Th (see Section 


6.1), thus keeping the ration H/h fixed. 


What is left now in order to turn the abstract estimate of Theorem 4.2 into a fully 
computable one is to specify s and v, given a DG solution ph. We will do so in the 
following paragraphs, finally using the knowledge that ph was computed using our DG 
discretization. 


Oswald interpolation. The form of the nonconformity indicator in Theorem 4.2 
already indicates how to choose s: it should be close to ph, in order to minimize fjnc, 
and it should be computable with reasonable effort. Both requirements are met by the 
Oswald interpolation operator, that goes back to j35| (in the context of a-posteriori error 
estimates; see also |26[ Section 2.5] and the references therein). Given any nonconforming 
approximation ph G Q^{Th) ^ we choose s G Rq(D) as a conforming reconstruc¬ 
tion of Ph by the Oswald interpolation operator Iqs : Q\{Th) G Rq(D) which 

is defined by prescribing its values on the Lagrange nodes u of the triangulation: we set 
Ios'[Ph]{v) '■= Ph{v) inside any t G and 

Ios\Ph]{v) '■= p\{v) for all inner nodes of Th and Ios\Ph]W) '■= 0 


ter- 


for all boundary nodes of t/j, where tJ) C denotes the set of all simplices of the fine 
triangulation which share as a node. 

We continue with the specification of v, which is a bit more involved. The only formal 
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requirement we have is for v to fulfill the local conservation property on each coarse 
element, although the diffusive flux estimator already gives a good hint on the specific 
form of V (namely, to be close to —A(/x)k£V/iP/i(/^)). A particular choice is given by the 
element-wise diffusive flux reconstruction (with respect to the fine triangulation) that was 
proposed in [26], which fulfills the local conservation property on the coarse elements if 
properly defined with respect to our DG discretization. 

Diffusive flux reconstruction. We will reconstruct a conforming diffusive flux ap¬ 
proximation Uh{^i) G iifdiv(f^) of the nonconforming discrete diffusive flux 0 

Hdivi^) in a conforming discrete subspace Vl{Th) C Hdiv{^), namely the Raviart- 
Thomas-Nedelec space of vector functions (see |26) and the references therein), defined 
for /c — 1 < / < A; by 

yki'^h) ■= G Fdiv(f^)|u|i G [Pz(t)]‘^-h®lPz(t) Vt G Th). 


See |261 Section 2.4] and the references therein for a detailed discussion of the role of 
the polynomial degree I, the properties of elements of Vl{Th) and the origin of the use 
of diffusive flux reconstructions in the context of error estimation in general. We define 


the diffusive flux reconstruction operator : V 
RhWh; At]], by locally specifying Ri[qh; fj] G Vl{Th). 

{Rhkh;tA-ne,q)L2^f, = KiQh, q; fj-) + bp{qh, g; tJ^) 


mrH) 

such that 


Ki^h)], fJ. ^ [qh ^ 


for all q G IPz(e) 


(14) 


and all e G and 


{Rh[qh;R],^q)j^ 2 ^i = -b\qh, q-, tJ.) - bl{q, qh] for all Vq G [IP/-i(t)]'^ (15) 




with q G IPz(t) for all t G r/j, where r? is given by the local discretization inside each coarse 
element and by ?? = 1 on all fine faces that lie on a coarse face. The next Lemma shows 
that this reconstruction of the diffusive flux is sensible for the DG solution as well as the 
reduced solution, since the reconstructions of both fulfill the requirements of Theorem 


Lemma 4.4 (Local conservativity). Let 1 G ^ ^ 

PhifJ-) G Q^{Th) be the DG solution of (|^ and Predip-) G QvediTn) the reduced so¬ 
lution of for a parameter pL G V and let Uhifx) := /i] G D^(r/i) and 

Ured(A*) := -R/i[A'red(A*); m] G V)(('rh) denote their respective diffusive flux reconstructions. 
It then holds that Uh{fx) and Ured{P') fulfill the local conservation property of Theorem 
4-31 i.e., 


(V-M/i(Ai),l)^2^-r = (/,l)i2^r= {'^■Ured{p-),L)p^ 2 ^T for all T gTh- 

Proof. We follow the ideas of |26| Lemma 2.1] while accounting for the coarse triangula¬ 
tion. Let G Q^(Tff) be an indicator for T G Th-, such that = 1 G and zero 
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everywhere else. It then holds with * = h, red, that 






= Ai) = (/,l)i2^r, 

for all T G Th^ where we used Green’s Theorem in the first equality, the definition of the 


diffusive flux reconstruction, (14), (15), and the definition of and in the second 
and the fact, that 1 G Qh'^ P* solves ([^ or respectively, in the third. 


□ 


Inserting the Oswald interpolation for s and the diffusive flux reconstruction for v in 
Theorem |4.2| then yields a locally computable energy estimate for the DG as well as the 
reduced solution. 

Corollary 4.5 (Locally computable energy norm estimate). Let p{^i) G Hq{Q,) be the 
weak solution of let PhifJ-) £ Q\0'h) DG solution of let p^g^{pL) G 

QiediTn) be the reduced solution of and let denote the diffusive flux reconstruction 
operator. Let the assumptions of Theorem (4.2) and Lemma (4.4) be fulfilled and let 
Ji, fi €V be arbitrary. It then holds, that 

WlpitJ-) - < ??(pft(At);/r,p, fi), 

|||p(/r) -pred(At)|||p: < ??(Pred(At);At,74, A)- 

with 






1/2 / 


and 


+ 


T&Th T£Th 

1/2 


\/a(M 


b)( E 


T&Th 


Vnci'1 • Vnci" •’^os[']i ! Pr (‘lA^) • Vi {Dhl', fx]), Vdfi'i P'') ' Vdli" i Dfi[', ffl, fl') ■ 

Local efficiency. The global efficiency of the abstract estimate was already shown in 
Lemma [4 .1 1 (again note, that 'yifx, JZ) = a{pi, Ji) = 1 if fx and /J coincide); see |251 Remarks 
4.2 and 4.3] for a discussion of the global efficiency of the abstract estimate. We also 
state a local efficiency of the local estimates using the knowledge of the discretization, the 
Oswald interpolation and the diffusive flux reconstruction. To do so we further localize 
our estimates with respect to the fine triangulation and apply the ideas of IHIESIEB]. 
We denote by < a proportionality relation between to quantities a and b in the sense 
that a <b :<;=> a < Cb, where the positive constant C depends on the space dimension, 
the polynomial degree k, the polynomial degree of /, the shape-regularity of Th and the 
DG parameters a and d, only. We additionally denote the set of all fine elements that 
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touch T G Th by := {t € Th\t CiT ^ 9}, the set of all fine faces that touch T by 
^ := {e G Jui I G : e H t / 0} and the weighted jump seminorms [[[•]]]. : P —?■ 
[H^{Th) -^R], ML f] and [[[-L,. - : V [H^{rh) -^R], 


for any subset T C Th-, all /x G P and q G H^{Th) by 

X 1/2 


fj.,T 


^ ||[[(A(/x)K£Vg) 
eSJ' 


Ue 


1^2,6/ 


and 


p,h>-T ’ 


'^bp{q,q;fi 

e^T 




1/2 


respectively. We analogously denote the set of all fine elements that touch t G by L 
and the set of all fine faces that touch t by and define 


cj := mine*, 
t&fT 


Cp := maxC* 
tern 


Ht ■= maxht, 


—T 

(jj 


iSr/ 


-( 


max oj, 


+ 2x1/2 




Kt := min hj, 


iSr/ 


Cj := maxC* 




C'L=maxf( max 

terj ^ * * 

for all T G 7h, where C* > 0 denotes the Poincare constant on t G T/i, defined analogously 
to Cp. 

Theorem 4.6 (Local efficiency of the locally computable energy norm estimate). With 
the notation and assumptions from Corollary \4.5\ let f be polynomial and max^g,-^ ht < 1. 
It then holds with * = h,i'ed, respectively, that 

'f'T /;;'r\l/2 


hlc (P* (m) ; /^) ^ {Cj/Ce) ^ . 




Cehr Mt^) - pM\\\jz,t 


+ Mf^) - pMlp,ji,j 

^df(p*(M);/^: A) ^ \/7 (m,A)\/7(m,A) lb(M) -p*(/^)IL,t 

+ MiJ-) - p*ifJ'Mp,-ii,Tf) 

for all coarse elements T G Th- 

Proof. We estimate each local estimator separately. It holds for the nonconformity esti- 


17 







mator, that 




N 1/2 


iSr/ 




T 


1/2 


t&Ti 

■yT /;-r\l/2r 


where we use the definition of i?ncb*(/^)] '^h equality and the arguments of 

|26[ Proof of Theorem 3.2] in the first and the definition of Cj, and in the second 
inequality. 

It holds for the residual estimator that 


: = (i) 


+ i|V- (A(/x)KeVp*(/x) +U*(/x))|| 


T P 


where we used the definition of ? 7 j’[u*(/x)] and the triangle inequality, which leaves us 
with two terms we will estimate separately. 

(i) The first term can be estimated as follows, using the definition of and the 
arguments of |24l Proposition 3.3] in the first and the definition of C'^ and and 
the fact that maxt^T^ hf < 1 in the second inequality: 


— 1/2 

(i) ;s E dv'iWf) - p.willed " 


<cjlv 'IIIp(m)-p.(m)IL,7. 


ter/ 


{a) The second term can be estimated as 


< 


[ii 


< 


^ Clhtci ut II [[(A(At)KeVp*(/x)) • 

t&Th e&rl 

+ Y - P*{P)lY,Ti ^ ' 

iSr^ 


rip. 


Il2 


using the definition of and the arguments of |26| Proof of Theorem 3.2] in the 
first and the definition of Cj, hx, Cj, and tT in the second inequality. 
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Applying the norm equivalence Q yields the desired result for the residual estimator. 
Finally, it holds for the diffusive flux estimator, that 


< \/7 (m,A)( ||(A(At)K£) ^^^(A(/i)k£Vp*(/ 2 ) +n*(/^)) 


t&ri 


2 y/2 

L2 J 


1/2 


< 


\/7(m,A) Qt 


■T1/2 


+ [b(M) 

- P* (m) III ;x,T + [[b(M) - P* (m)]]]p,m, ) > 


where we used the definition of 7jf[p* (/.*)] and and the parameter equivalence from 
(|^ in the first, |25[ Lemma 4.12] in the second and the definition of cj and in the 
third inequality. Applying the norm equivalence Q again finally yields the desired result 
for the diffusive flux estimator. □ 


5 On-line enrichment 

As mentioned in Section [3] there are two drawbacks of classical RB methods in the con¬ 
text of parametric multi-scale problems, stemming from the fact that dimQ^ roughly 
scales with 0{e~^): (i) expensive high-dimensional evaluations of global quantities dur¬ 
ing reduction and orthogonalization and (ii) expensive high-dimensional inversions of 
(|^ during the basis generation. The first shortcoming was addressed by the LRBMS 
method introduced in |36[ [8] and finalized by the local error estimator introduced in the 
previous section, that can be off-line/on-line decomposed with a computational complex¬ 
ity depending only linearly on |r/i| (not shown) while the original estimator presented 
in [8] required the computation of |7A| global Riesz-representatives for each snapshot. 
But also the LRBMS method suffers from the second shortcoming, namely that the basis 
generation (for instance using an adaptive Greedy procedure) requires the global high¬ 
dimensional problem ([^ to be solved for a possibly large number of parameters out of a 
finite training set 'Ptrain 7 V. In the context of parametric multi-scale problems, however, 
solving ^ may be prohibitively expensive and one may only have the resources to do so 
for a very limited amount of parameters, if at all. Such an RB space constructed out of 
only very few solution snapshots is usually insufficient for nearly all parameters fx € V 
that were not included in the basis generation. 

To address this shortcoming we relax the notion of a strict off-line/on-line splitting of 
the computation in the classical sense. While the computational complexity of the on-line 
phase must not depend on any global high-dimensional quantities (namely dimQ^) for 
RB methods, we allow for high-dimensional but local computations (of order dimQ^’^, 
with T G 7h) hr the context of the LRBMS method. The idea is as follows: during the 
off-line phase we initialize the local reduced bases with a DG basis of order up to 
kn G N with respect to the coarse elements T G Th-, thus ensuring that any reduced 
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Algorithm 5.1 Discrete weak Greedy algorithm used in the LRBMS method 

Input: ONB, kn £ N, 'Ptrain C P, Agreedy > 0, Ngreedy £ N 
Output: A local reduced basis for each coarse element T £ 7h- 
Initialize the local reduced bases with the coarse DG basis: 

for all T € 7h do 

^t(o) ^ onb({DG shape functions of order up to ku w.r.t. T}) 

end for 

n ■«— 0 

while max rj{pred{fj-)-, 11,71, II) > Agreedy, withpred(p) G Qred(TH)^"^ Solviug ® 
n-ePtrain ' 

and n < Ag^eedy do 

Compute all reduced quantities w.r.t Qred(Th)^"^- 

Find the worst approximated parameter, with p^edili) £ Qred(Tr/)^"^ solving §: 

Mmax ^ argmax ri{p,edifJ.)-, fJ., /r, A) 


© 


span(^'^™) 


TGTh 


nS'Pti, 

Extend the local reduced bases, with ph (p) £ Qh solving 

for all T £ Th do 

end for 

gred(rH)'"+^^ 
n n + 1 

end while 

iTeTn 


© 


TSTff 


Ph(/x)|y}) 

span(^ 


return [HP 


solution is at least as good as a DG solution on the coarse triangulation. We then carry 
out a discrete weak Greedy algorithm based on the error estimator dehned in Gorollary 
4.5 while allowing only for a limited amount of global solution snapshots, Aynax £ N, and 
extend the local bases with these snapshots using an orthonormalization algorithm ONB 
locally on each T G Th- This procedure is summarized in Algorithm |5.1| 

During the on-line phase, given any € V, we compute a reduced solution pi.ed(/^) £ 
QrediTu) and efficiently assess its quality using the error estimator. If the estimated 
error is above a prescribed tolerance, Aoniine > 0, we start an intermediate local enrich¬ 
ment phase to enrich the reduced bases in the SEMR (solve —)• estimate —)• mark —)• 
rehne) spirit of adaptive mesh rehnement (the procedure is summarized in Algorithm 
5.2): we first compute local error indicators rj^fi) for all T G Th, such 


that r]{-; /x, fi, fiY < T.T£Th ^ ’ dehned as 






V7(AbA0 Anc(-; fj-f + vl (•; fj-? + 




Adf(-;At,A)^ (16) 


and mark coarse elements Th T Th for enrichment, given a marking strategy MARK. For 
each marked T G Th we solve 


blYpl'{t^):Qh;t^) = IhiQh) 


— Tsi 


for all qh G 


(17) 


on an overlapping domain Ts D T with the insufficient reduced solution pred(M) ns dirich- 
let boundary values on dTs to obtain an updated detailed solution p^'^(/x) G 
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Here ^ extensions of Q^i'^'h)^ i respectively, to the over- 

sampled domain T 5 while additionally encoding Pred(/^) as dirichlet boundary values. 
We then extend each marked local reduced basis by performing an orthonormalization 
procedure on with respect to the existing local basis and update all reduced 

quantities with respect to the newly added basis vector. We finally compute an updated 
reduced solution in the updated coarse reduced space and estimate the error again. We 
repeat this procedure until the estimated error falls below the prescribed tolerance Aoniine 
or until the prescribed maximum number of iterations, A^online £ N, is reached. Possible 
choices for ONB and MARK are given in Section 6.2 Note that first steps in the direction 
of on-line enrichment were published in [9]. 


Algorithm 5.2 Adaptive basis enrichment in the intermediate local enrichment phase 

Input: MARK, ONB, Pved(fJ-), P., Aoniine > 0, -^online ^ ^ 

Output; Updated reduced solution and local reduced bases. 

^t{0 ) ^ yy ^ 

n <— 0 

while J7(pred(/.t); A) > Aoniine and n < Aoniine do 

for all T G Th do 

Compute local error indicator {py.ed{lJ-)\ fa, Ji, fi) according to ( |16[ ). 
end for 
Th ^ MARK(r//) 
for all T G Th do 

Solve forpl^in) G QUt')- 

end for 

AedlTh) 

Update all reduced quantities w.r.t Qred(TH)^"^^A 
Solve 0 for Prediu) ^ Qred(Tr/)^"^^^ 
n •<— n -I- 1 

end while 


^7^(n+l)\ 


Qred(TH)^"+^^ ^ 0T67>f Spanspan(<A^'"^) 


return p^ediu), } 


TGTh 


Remark 5.1. It is also possible to use other methods to compute (or approximate) Phift) 
during the Greedy procedure, in particular other domain decomposition or multi-scale 
methods. If the resulting approximation does not fulfill the local conservation property of 


Lemma 4-4 • however, the estimator would have to he replaced during the Greedy algorithm, 
for instance in the spirit of m- But during the on-line phase the estimator would he 
valid for any reduced solution, as long as the basis is initialized with at least a constant 
function. In particular one could use variants of the MsFEM (see ]22f). the HMM (see 
m) or in particular the DG-HMM (see mw during the Greedy procedure to generate a 
coarse reduced basis with approximation properties of order H. Fine scale features of the 
solution would then be adaptively recovered during the on-line enrichment phase, if and 
where needed. 


Remark 5.2. The computation of the local error indicators in Algorithm \5.2\ can be effi¬ 
ciently off-line/on-line decomposed (not shown here). Once a set of subdomains has been 
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marked, the enriehment ean be earried out in parallel without any eommunieation. For 
the update of the reduced quantities only local information and the information on one 
layer of neighboring fine grid cells is needed. The on-line phase, however, requires infor¬ 
mation on Th (in particular the number of coarse elements and neighboring information), 
in contrast to traditional RB methods. 

Remark 5.3. Our choice of the Greedy algorithm and the adaptive on-line enrich¬ 
ment covers a wide range of scenarios. Disabling the on-line enrichment (by setting 
Nonline = Oj and choosing any suitable Ag^eedy o.'n-d Ngreedy yields the standard Greedy 
basis generation, well known in the RB context. Setting Ngreedy = 0 and kn = 1, on the 
other hand, disables the Greedy procedure and merely initializes the reduced bases with 
the coarse DG basis of order one. This is of particular interest in situations where the 
computation of solutions of the detailed problem during the Greedy procedure might be too 
costly. In that setup nearly all work is done in the adaptive on-line enrichment phase. 
Many other variants are possible, e.g. other local boundary conditions, several mark¬ 
ing strategies MARK or orthonormalization algorithms ONE or other stopping criteria; one 
could also limit the number of intermediate snapshots added to the local bases. Depending 
on these choices the resulting method is then close to existing DD methods (i.e., a DD 
method with overlapping subdomains, see m or multi-scale methods (i.e., the adaptive 
iterative multi-scale finite volume method mi)- 

The LRBMS method with the proposed adaptive on-line enrichment strategy is now 
snitable for a far wider range of circnmstances than standard RB methods or the pre- 
vionsly pnblished variant of the LRBMS method. As mentioned before it can now be 
applied if the compntational power available for the off-line phase is limited by time- or 
resonrce constraints. It can also be applied if the set of training parameters given to the 
Greedy Algorithm was insnfhciently chosen or even if on-line a solntion to a parameter is 
reqnested that is ontside of the original bonnds of the parameter space. In general, the 
on-line adaptive LRBMS method can be applied whenever the basis that was generated 
dnring the off-line phase tnrns ont to not be snfhcient for what is reqnired dnring the 
on-line phase. 

6 Numerical experiments 

In this section we investigate the performance of the error estimator in the context 
of the DG discretization defined in Section 12.21 as well as in the context of the the 
LRBMS method defined in Section We also investigate the performance of the on-line 
enrichment procednre we proposed in the previons section. 

We nsed several software packages for the implementation: everything concerning the 
discretization was implemented within the high performance G++ software framework 
DUNE da [TT] while everything related to model rednction was implemented based on 
the pyMOR package m- Data fnnctions, container and linear solvers were implemented 
within dune-stuff im, discrete function spaces (based on the discretization modules 
dune-fern |19| and dune-pdelab [T]), operators and products for the discretization as 
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well as the error estimator were implemented within dune-gdt |49| . The coarse triangu¬ 
lation was implemented in dune-grid-multiscale |50| while the Python bindings form¬ 
ing the bridge between DUNE and pyMOR were implemented in dune-pymor j48j . Finally a 
high-level solver for linear elliptic (and possibly parametric) PDFs was implemented in 
dune-hdd |51) and exposed to pyMOR, where everything related to model reduction was 
implemented. 

The fine triangulations Th used are conforming refinements of triangular grids, rep¬ 
resented by instances of ALUGrid< 2, 2, simplex, conforming > (see |18)). In the 
following experiments we choose a DG space {t'^) locally on all coarse elements 
T G Th- The resulting discretization thus coincides with the one defined in |27[ [26] 
(though in general other choices are possible). All coarse triangulations Th used consist 
of squared elements (though arbitrary shapes are possible). 

Most figures in this publication were created or arranged using Ti^Z |55[ l54] and 
pgfplots [28], the colorblind safe colors in Figures]^ and were selected with Color Brewer 

HSj. 


6.1 Study of the a-posteriori error estimator 

To study the convergence properties of our estimator we consider two experiments. The 
first one serves as an academic example and as a comparison to the work of |25[ |26] . The 
second experiment demonstrates the efficiency of the estimator in realistic circumstances. 
In both experiments we compute estimator components rjf^ := YIt&Th 6*“^ ^ 
r, df and the estimator r] as defined in Corollary |4.5| using a 0th order diffusive ffux recon¬ 
struction (Z = 0 in (14), (|15|)). If no analytical solution p(/i) is available we approximate 
the discretization error |||p(/x) — p{tJ-)\\\-p by substituting p(/i) for a discrete solution on 
a finer grid and by computing all integrals using a high order quadrature on the finer 
grid. We denote the efficiency of the estimator, T]{ph{fi)-, /x,/!, jj,)/\\\p{fj,) — PhifJ')\\\jj; > 1; 
by “eff.” and the average (over all refinement steps) experimental order of convergence 
of a quantity by “order”. 

Academic example. We consider ([^ on 12 = [—1,1]^ with a parameter space V = 
[0.1,1], Ke = id, f{x, y) = cos(|7rx) cos(^7r?/), A(x, y; /x) = 1-|-(1—/x) cos(|7rx) cos(^7ry) 
and homogeneous Dirichlet boundary values. We study the components of the estima¬ 
tor as well as its efficiency in several circumstances, i.e., for different parameters /x, /x, 
fi € V and triangulations Th and Th- We choose Th just as in |25| Section 8.1] and 
begin with pi = Ji = fd = 1, thus reproducing the nonparametric example studied in 


WA 

II|P(/X) -Ph(/X)|||j;j- 

ync(-;/x) 


ydf(-;/x, A) 

eff. 

128 

3.28-10“^ 

1.66-10"^ 

5.79-10-" 

3.55-10"^ 

3.36 

512 

1.60-10“^ 

7.89-10"^ 

2.90-10"^ 

1.76-10"^ 

3.40 

2,048 

7.78-10-2 

3.91-10“^ 

1.45-10"^ 

8.73-10"^ 

3.49 

8,192 

3.47'10“^ 

1.95-10"^ 

7.27'10"^ 

4.35-10"^ 

3.91 

order 

1.08 

1.03 

1.00 

1.01 

- 


Table 1: Discretization error, estimator components and efficiency of the error estimator 
for the academic example in Section 6.1 with \Th\ = 1 and /x = /x = /x = 1. 
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A = 1 

A = 0.1 

Iml / [ThI 


y{-',y-T, A) 

eff. 

A) 

n{-\y-,V; A) 

eff. 

128/ 2x2 

2.89-10“^ 

8.10'10"" 

2.47 

3.16-10"" 

7.71-10"" 

2.35 

512/ 4x4 

7.26-10“^ 

3.27-10"^ 

2.04 

1.56-10"^ 

3.08-10"^ 

1.92 

2,048/ 8 x8 

1.82-10“^ 

1.45-10"^ 

1.86 

7.74-10"^ 

1.35-10"^ 

1.73 

8,192/16 X 16 

4.54-10“® 

6.76-10"^ 

1.95 

3.85-10"^ 

6.26-10"^ 

1.80 

order 

2.00 

1.20 

- 

1.01 

1.21 

- 


Table 2: Selected estimator components, estimated error and efficiency of the error esti¬ 


mator for the academic example in Section 6.1 with and Th simultaneously 
rehned, /.r = 1 and two choices of fi. Note that the estimator components r/nc 
and ? 7 df are not influenced by Th-, the estimator components r/nc and r/r are not 
influenced by fi and the discretization error is not influenced by either. Thus 
only 7 /r, 7] and its efficiency are given for /i = 1 and only r/df, r] and its efficiency 
are given for /i = 0.1 (the other quantities coincide with the ones in Table l|. 


|r.| / \Th\ 

Mfj-) -PhT)\\\^ 

Vnci-T) 

A) 

eff. 

128 / 2x2 

3.81-10~" 

1.82-10~" 

1.18-10“ 

3.10 

512/ 4x4 

1.87-10-^ 

8.57-10-2 

5.00-10-1 

2.67 

2,048 / 8 X 8 

9.08-10-2 

4.22-10"^ 

2.29-10“! 

2.52 

8,192/16 X 16 

4.05-10-2 

2.11-10-2 

1.10-10“! 

2.71 

order 

1.08 

1.03 

1.14 

- 


Table 3: Discretization error, selected estimator component, estimated error and effi¬ 


ciency of the error estimator for the academic example in Section 6.1 with 


and Th simultaneously rehned, 1 and = fi = 0.1. Note that rj,- and rj^i 
are not inhuenced by JI and coincide with Table 


|251 Section 8.1] (since A = 1 and all constants involving a and 7 are equal to 1). For 
this specihc choice of parameters an exact solution is available (see |25l Section 8.1]). 
In this conhguration, and coincide with their respective counterparts dehned in 
(251126] while is directly inhuenced by the choice of the coarse triangulation and the 
parametric nature of A (entering cj). Choosing Th = ^ (the coarse grid conhguration 
with the worst efficiency), we observe results similar to |25l Table 1] for r/df and r^nc in 
TableIn contrast, r/r shows only linear convergence while the residual estimator in |25| 
Table 1] converges with second order. Overall, the efficiency of the estimator ry is around 
3.5 (for hxed \Th\ = 1) while the efficiency of the estimator in |25[ Table 1 ] is around 
1 . 2 . We can recover the superconvergence of r/r, however, by rehning Th along with Th 
(thus keeping the ratio H/h hxed), see the left columns of Table In order to make 
the estimator off-line/on-line decomposable in the parametric setting, jl has to be hxed 
throughout the experiment. Choosing fi = 0.1 has no negative impact on the efficiency 
of the estimator, as we observe in the right columns of Table It is often desirable to 
additionally hx the error norm throughout the experiments. Choosing /J = 0.1 we still 
observe a very reasonable efficiency in Table 

Multi-scale example. We consider (§ on D = [0, 5] X [0,1] with f{x, y) = 2 ■ 10^ if 
{x,y) G [0.95,1.10] x [0.30,0.45], f{x,y) = -1 • 10 ^ if {x,y) G [3.00,3.15] x [0.75,0.90] or 
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|rd / \Th\ 

IIIp(m) -Ph(M)lllir 

»7nc(-;M) 

yr(-;M) 

Vdi{-;fx,P) 

eff. 

16,000 / 25 X 5 

7.49-10~" 

2.13-10*^ 

1.88-10"° 

9.66-10"" 

4.14 

64,000/ 50 x 10 

4.52-10~^ 

1.46-10° 

7.05-10"^° 

6.05-10"" 

4.58 

256,000/100 X 20 

2.58-10~^ 

1.02-10° 

7.44-10"“ 

3.85-10"" 

5.44 

1,024,000/200 X 40 

1.26-10“^ 

7.20-10"^ 

2.00-10"^° 

2.49-10"" 

7.70 

order 

0.86 

0.52 

1.07 

0.65 

- 


Table 4: Discretization error, estimator components and efficiency of the error estimator 


for the mnlti-scale example in Section 6.1 with Th and Th simnltaneonsly refined 
and fx = JZ = p. = 1. Note that rj,- shonld be close to zero (since / is piecewise 
constant, compare Fignre Q and snffers from nnmerical inaccnracies (ignoring 
the last refinement wonld yield an average order of 2.33 for r/r). 


(x, y) G [4.25,4.40] x [0.25,0.40] and 0 everywhere else, A(x, y; /x) = 1 + (1 — /x)Ac(x, y), 
Ate = ^id, homogeneons Dirichlet bonndary valnes and a parameter space V = [0.1,1]. 
On each t € Th, is the corresponding 0th entry of the permeability tensor nsed in the 
first model of the 10th SPE Comparative Solntion Project (which is given by 100 x 20 
constant tensors, see |5^) and Ac models a channel, as depicted in Fignre The right 
hand side / models a strong sonrce in the middle left of the domain and two sinks in 
the top and right middle of the domain, as is visible in the strnctnre of the solntions 
(see Fignre [^. The role of the parameter fi is to toggle the existence of the channel Ac. 
Thns X{h)k rednces to the above mentioned permeability tensor for /x = 1 while = 0.1 
models the removal of a large condnctivity region near the center of the domain (see 
the first row in Fignre]^. This missing channel has a visible impact on the strnctnre 
of the pressnre distribntion as well as the reconstrncted velocities, as we observe in the 
left colnmn of Fignre]^ With a contrast of 10® in the diffnsion tensor and an e of abont 
|D|/2,000 this setnp is a challenging heterogeneons mnlti-scale problem. 

While we observe in Table that the convergence rates of the estimator components 
are not as good as in the experiment stndied before, the estimator shows an average 
efficiency of 5.5, which is quite remarkable considering the contrast of the data functions. 
Fixing pL = p = 0.1 (to obtain a fully off-line/on-line decomposable configuration with 
a fixed error norm) yields an average efficiency of 10.2 which is still very reasonable. In 
addition we observe a good agreement between the spatial distribution of the error and 
the estimator indicators in Figure]^ 



Figure 1: Data functions of the multi-scale example in Section |6.l| on a triangulation with 
|t/i| = 16,000 elements: location of the channel function Ac (left) and plot of 
the force / (right) modeling one source (black: 2-10®) and two sinks (dark gray: 
—1-10^, zero elsewhere). 
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H = 1 


n = 0.1 



Figure 2: Data functions and sample solutions of the multi-scale example in Section 6.1 on 
a triangulation with |t/i| = 16,000 elements for parameters /i = 1 (left column) 
and /X = 0.1 (right column). In each row both plots share the same color 
map (middle) with different ranges per row. From top to bottom: logarithmic 
plot of A(/x)k (dark: 1.41-10“^, light: 1.41-10^), plot of the pressure Ph(/^) 
(solution of Q, dark: —3.92-10“^, light: 7.61-10—1, isolines at 10%, 20%, 
45%, 75% and 95%) and plot of the magnitude of the reconstructed diffusive 
flux fj] (defined in (14) and (15), dark: 3.10-10“®, light: 3.01-10^). 


Note the presence of high-conductivity channels in the permeability (top left, 
light regions) throughout large parts of the domain. The parameter dependency 
models a removal of one such channel in the middle right of the domain (top 
right), well visible in the reconstructed Darcy velocity fields (bottom). 


n = 1.0 


n = 0.1 



Figure 3: Spatial distribution of the relative error contribution (top), 
\\\p{tJ-) - PhitJ-)\\\-p^T/\\\p(t^) ~ Ph{f^)\\\-p:^ and the relative estimated error 

contribution (bottom), rf {ph{p)] P-)/{Y.T£Th (Phif^): 

for all T G Th , for the multi-scale example in Section |6.1| with | Th | = 25 x 5 
and Jl = fi = 0.1 for parameters = 1 (left column, light: 2.26-10“^, dark: 
3.78-10“^) and pi = 0.1 (right column, light: 4.02-10“^, dark: 3.73-10“^). 
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6.2 Adaptive on-line enriehment 


To demonstrate the proposed adaptive on-line enrichment Algorithm |5.2| and the flexi¬ 
bility of the LRBMS we study two distinct circumstances. We first consider a smooth 
academic example where we disable the greedy procedure and build the reduced bases 
only by local enrichment. The second example is again a multi-scale one with global 
channels in the permeability and we allow for very few global solution snapshots and 
adaptively enrich afterwards. 


For the orthonormalization algorithm ONB in the Greedy Algorithm 5.1 we use the 
stabilized Gram Schmidt procedure implemented in pyMOI0 with respect to the scalar 
product given by {p, q) i—)• b'^{p, q; Ji) YleeTT bp {p, q] Jl) on each T G Th- In contrast to 
other possible basis extension algorithms (e.g. using a proper orthogonal decomposition) 
the Gram Schmidt basis extension yields hierarchical local reduced bases. This is of 
particular interest in the context of on-line enrichment, since we do not need to update 
reduced quantities w.r.t existing basis vectors after enrichment. We always initialize the 
local reduced bases with the coarse DG basis of order up to one by setting kn = 1 in the 
Greedy Algorithm |5.1| We choose the same orthonormalization algorithm in the adaptive 
basis enrichment Algorithm |5.2| and use several marking strategies for MARK, depending 
on the circumstances (detailed below). Regarding the overlap for the local enrichment 
we always choose the overlapping subdomains D T to include T and all subdomains 
that touch it, thus choosing an overlap of 0(H) as motivated in |32) . 

Since the error of any reduced solution |||p(/.t) — Pred(/^)|||p; can not be lower than the 
error of the corresponding detailed solution |||p(/.t) — A'/i(M)|||p: we always choose Aoniine in 
Algorithm |5.2| to be slightly larger than max^gp^^j.^^ fi, JZ, jl) in our experiments 

(see below), where 'Ponline C "P is the set of all parameters we consider during the on-line 
phase. This is only necessary since we do not allow for an adaptation of t/j; combining 
our on-line adaptive LRBMS with the ideas of |58| would overcome this restriction. 

Academic example. We again consider the academic example detailed in Subsection 
6.1 on fixed triangulations with |r/i| = 131,072 and \Th\ = 8 x 8 and choose the set of 


on-line parameters Ponline to consist of 10 randomly chosen parameters pg,... , 7 X 9 G P. 
We set A^greedy = 0 and kn = 1, thus disabling any Greedy search and initializing the 
local bases with the coarse DG basis of order up to one (consisting of 4 shape functions). 
In this setup max^gp^^j.^^ ri{ph{pi)-, p,/I, ji) = 2.79 • 10“^ and we choose Aoniine = 5 • 10“^ 
in Algorithm |5.2[ 

We begin by choosing MARK such that Th = Th (ah coarse elements are marked; 
denoted by uniform in the following). For each parameter p G Ponline this results in 
a method that is similar to DD methods with overlapping subdomains. In contrast to 
traditional DD methods, however, we start with an initial coarse basis and perform a 
reduced solve before each iteration which helps to spatially spread information. As we 
observe in Figure]^ top left, it takes four enrichment steps to lower the estimated error 
for the first on-line parameter pg below the desired tolerance and another two enrichment 
steps for the next on-line parameter pui (which is max^g-p^^^jj^^ p). With uniform marking 


Tttp: //docs .pymor. org/en/0.2 . x/_modules/pymor/la/grain_schinidt .html#gram_schinidt 
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uniform 


doerfler_age(1/3,4) 


10 " 


<a, 

li 


Sr- 10“ 



dimQred(7H) 

uniform_doerfler_age(10,1/3, 4) 


dim QrodCTf/) 



^online 


A fig = 0.43708. 
-m-lii = 0.95564. 
- /i2 = 0-75879. 
A #^3 = 0.63879. 
0 M4 = 0.24041. 
ic fig = 0.24039. 
-B-Me = 0.15227. 
X M7 = 0.87955 . 
+ /Tg = 0.64100. 
O Mg = 0.73726. 


Figure 4: Estimated error evolution during the adaptive on-line phase for the academic 
example in Section 6.2 with \Th\ = 64, A^gj-eedy = 0, kn = 1, Aoniine = 5-10“^ 
(dotted line) and Ji = p, = 0.1 for several marking strategies: uniform marking 
of all subdomains (top left), combined Dorfler marking with 0doerf =1/3 and 
age-based marking with Aage = 4 (top right) and additional uniform marking 
while r/(pi.ed(/-i); A) > ^uniAoniine with ^ nni = 10 (bottom left). With each 
strategy the local rednced bases are enriched according to Algorithm |5.2| for 
each snbsequently processed on-line parameter /Xq, ..., /Xg (bottom right). 
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this increases the local basis sizes from four to ten on each coarse element. The resulting 
coarse reduced space of dimension 640 is then sufficient to solve for the next four on-line 
parameters / 22 ,..., /X 5 without enrichment. One additional enrichment phase is needed 
for /Xg (which is min^g-p^^j.^^ /x) and none for the remaining three on-line parameters. Note 
that while this uniform marking strategy may be optimal in the number of enrichment 
steps it takes to reach the desired error for all on-line parameters it also leads to an 
unnecessarily high-dimensional coarse reduced space (of dimQred(TH') = 704) and a high 
work-load in each enrichment step. 

Another popular choice in the context of adaptive mesh refinement is a Dorffer marking 
strategy (see |42) and the references therein), where we collect those coarse elements that 
contribute most to ^doerfX^TsTff ^ Th, for a given user-dependent 

parameter 0 < ^doerf < 1 - In addition, similar to na |29], we count how often each 
T G Th was not marked and mark those elements the “age” of which is larger than 
a prescribed A'age £ N (resetting the age count of each selected element). We denote 
this marking strategy by doerf ler_age(0doerf)-^age) ■ We found that a combination of 
^doerf =1/3 and A^age = 4 yielded the smallest overall basis size (of 572), compared to 
other combinations of 0doerf and A"age and the uniform marking strategy. The number 
of elements marked per step range between five and 52 (over all on-line parameters and 
all enrichment steps; 23 steps in total) with a mean of 14 and a median of ten. Of these 
marked elements between one and 44 have been marked due to their age in 12 of these 23 
steps (with an average of 12 and a median of eight, taken over only those 12 steps where 
elements have been marked due to their age). We observe in Figure]^ top right, that 
the general behavior of the method with this marking strategy is similar to the one with 
uniform marking, with some commonalities and differences worth noting. First of all it 
also takes three enrichment phases to reach the prescribed error tolerance, and for the 
same parameters /Xq, and /Xg as above. But each of these enrichment phases naturally 
need more steps and large improvements can usually be observed after a lot of elements 
have been marked due to their age count (see for instance the fifth enrichment step for 
/Xq or Hi). In addition we observe that the estimated error for a parameter sometimes 
increases, in particular in the very beginning (see the first four steps for /Xq, the fourth 
step for Hi or th® first step for /Xg). This is not troublesome since we can only expect 
a strict improvement in the energy norm induced by the bilinear form that is used in 
the enrichment. This shows, never the less, that there is still room for improvement, 
although using the doerfler_age marking we could reach a significantly lower overall 
basis size than using the uniform marking (572 vs. 704). 

We propose a combination of the two strategies, namely a uniform marking while 
the estimated error is far away from the desired tolerance, i.e., r 7 (pi.ed(M)j A) A) > 
0uni ^online for some 0uni > 0, followed by a Dorfler and age-based marking as detailed 
above. We denote this marking strategy by uniform_doerfler_age(0uni;^doerf; IWge)- 
As we observe in Figure]^ bottom left, this marking strategy combines advantages of both 
previous approaches, recovering the rapid error decrease of the uniform marking strategy 
far away from the desired tolerance (see the first step for /Xq and Hi) while yielding the 
smallest overall basis size of 530 (using a factor of 0uni = 10) due to the doerfler_age 
marking strategy. The smoothness and symmetry of the problem is reflected in the 
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Figure 5: Spatial distribution of the final sizes of the local reduced 
bases, (light: 7, dark: 11), for all T G Th after 

the adaptive on-line phase for the academic example in 
with Q = [—1,1]^, \Th\ =8 x8 and the 


Section 


6.2 


uniforin_doerfler_age(10,1/3,4) marking strategy (see 
Figure 1^ bottom left). 


spatial distribution of the final local basis sizes (see Figure]^. 

Multi-scale example. We again consider the multi-scale example detailed in §6.1| on 
fixed triangulations with |t/i| = 1,014,000 and \Th\ = 25 x 5 and choose the set of on-line 
parameters 'Ponline to consist of the same 10 randomly chosen parameters /Xq, ..., pg G "P 
as in the previous example. In this setup max^gp^^jj^^ ? 7 (p/i(/x); p, p,/t) = 2.66 and we 


choose Aoniine = 2 in Algorithm 5.2 


We first set Agreedy = 0 and kn = 1 (thus disabling any Greedy search in the off-line 
phase and initializing the local bases with the coarse DG basis of order up to one); in 
the on-line phase we use the uniform marking strategy (see above). As we observe in 
Figure]^ top, it takes 129 enrichment steps to lower the estimated error below the desired 
tolerance for the first on-line parameter /Tq. This is not surprising since the data functions 
exhibit strong multi-scale features and non-local high-conductivity channels connecting 
domain boundaries (see Figure]^. After this extensive enrichment it takes 12 steps for 
Hi and none or one enrichment steps to reach the desired tolerance for the other on-line 
parameters. The resulting coarse reduced space is of size 10,749 (with an average of 86 
basis functions per subdomain), which is clearly not optimal. Although each subdomain 
was marked for enrichment, the sizes of the final local reduced bases may differ since the 
local Gram Schmidt basis extension may reject updates (if the added basis function is 
locally not linearly independent). As we observe in Figure]^ left, this is indeed the case 
with local basis sizes ranging between 24 and 148. 

To remedy the situation we allow for two global snapshots during the off-line phase 
(setting Agreedy = 2 , Ptrain = { 0 . 1 , 1 }) and use the adaptive unif orm_doerfler_age 
marking strategy (see above) in the on-line phase. With two global solution snapshots 
incorporated in the basis the situation improves significantly, as we observe in Figure 
bottom left, and there is no qualitative difference of the evolution of the estimated error 
during the adaptive on-line phase between the academic example studied above and this 
highly heterogeneous multi-scale example (compare Figure]^ bottom left). In total we 
observe only two enrichment steps with uniform marking (see the first two step for Ho)- 
The number of elements marked range between 11 and 110 (over all on-line parameters 
and all but the first two enrichment steps) with a mean of 29 and a median of 22 . Of 
these marked elements only once have 87 out of 110 elements been marked due to their 
age (see the last step for Hi)- Overall we could reach a significantly lower overall basis 
size than in the previous setup (1,375 vs. 10,749) and the sizes of the final local bases 
range between only nine and 2 (compared to 24 to 148 above). We also observe in Figure 
right, that the spatial distribution of the basis sizes follows the spatial structure of 
the data functions (compare Figures mi- which nicely shows the localization qualities 
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Figure 6: Estimated error evolution during the adaptive on-line phase for the multi¬ 


scale example in ^6.2 with \Th\ = 125, kn = 1, Aoniine = 2 (dotted line), 
/I = /I = 0.1 for different on-line and off-line strategies: no global snap¬ 
shot (Greedy search disabled, A^greedy = 0) during the off-line phase, uniform 
marking during the on-line phase (top) and two global snapshots (Greedy 
search on T^train = {0.1,1}, iVgreedy = 2) and combined uniform marking 
while r/(pred(/^);/^) A) > ^uni^^online with 0uni = 10, Dorfler marking with 
^doerf = 0.85 and age-based marking with Aage = 4 (bottom left); note the 
different scales. With each strategy the local reduced bases are enriched ac¬ 
cording to Algorithm |5.2| while snbseqnently processing the on-line parameters 
fj,Q,..., fj,g (bottom right). 
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Figure 7: Spatial distribution of the final sizes of the local reduced bases, for all 
T G 7//, after the adaptive online phase for the multi-scale example in §6.2| 
with n = [0,5] X [0,1], \Th\ = 25 X 5 for the two strategies shown in Figure 
1^ no global snapshot with uniform enrichment (left, light: 24, dark: 148) 
and two global snapshots with adaptive enrichment (right, light: 9, dark: 20). 
Note the pronounced structure (right) reflecting the spatial structure of the 
data functions (compare Figures and |^. 


of out error estimator. 


7 Conclusion 


In this contribution we equipped the localized Reduced Basis multi-scale method with 
a rigorous and efficient localized a-posteriori error estimator and a new adaptive on-line 
enrichment procedure. The LRBMS method was originally introduced in |36l [8] with a 
standard residual based a-posteriori error estimator for the reduced solution against the 
high-dimensional solution that could not give local error information and was computa¬ 
tionally very costly during the off-line phase of the computation. In addition the method 
was restricted to the SWIPDG discretization. We extended the original LRBMS method 
to locally allow for arbitrary discretizations of at least first order within each subdomain. 
In addition, we proposed a new error estimator for the high-dimensional as well as the 
reduced solution against the weak solution in the spirit of |25l |26] , based on a local con¬ 
servative flux reconstruction. The proposed estimator yields a guaranteed upper bound 
(involving no unknown constants), it is locally efficient and can be computed using local 
information only. The estimator is efficiently off-line/on-line decomposable and allows to 
estimate the error in parameter dependent norms. 

Using the local information of the estimator we proposed a new adaptive on-line en¬ 
richment strategy, where we extend the local reduced bases by solutions of local corrector 
problems posed on overlapping subdomains. This is done during the on-line phase, thus 
deviating from the strict off-line/on-line separation of traditional RB methods, but only 
on those subdomains that have been selected by the estimator. This strategy allows to 
guarantee the quality of the reduced solution during the on-line phase, even if an insuf¬ 
ficient reduced basis had been prepared in the off-line phase (in contrast to traditional 
RB methods where the reduced basis is fixed for the whole on-line phase). 

We provide numerical experiments to demonstrate the performance of the proposed 
estimator in the parametric setting for the high-dimensional discretization for an aca¬ 
demic as well as a highly heterogeneous multi-scale example, where the estimator proves 


to be very efficient (see Section 6.1). We also demonstrate the performance of the newly 


proposed adaptive on-line enrichment strategy for both examples. It is particular note- 
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worthy that only very few global solution snapshots were needed during the off-line phase 
for the multi-scale example to sufficiently prepare a reduced basis that was then adap¬ 


tively enriched during the off-line phase (see Section 6.2). For the academic example 


of smooth data functions without any multi-scale features no global solution snapshots 
were required at all. 

First findings regarding the new estimator have been published in |45| while first steps 
in the direction of on-line enrichment have been discussed in j9]. 
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